», : 
Se ea 3 
. m= ate 
a “a 7 —, f : 
” ay, : : 
b e 
. ae < 3 
* aa : & 
> -* Zs 
, 
2 





Th, liaee ie 


SOLUTION OF FULL POTENTIAL EQUATION FOR THREE 
_ DIMENSIONAL SUBSONIC FLOW: BY A 
FIELD INTEGRAL METHOD 


by 
KUMAR ALUK | 





‘K, 98s | 


‘DEPARTMENT OF AERONAUTICAL ENGINEERING | 


— INDIAN INSTITUTE OF TECHNOLOGY, KANPUR 


JANUARY, 1989 


SOLUTION OF FULL POTENTIAL EQUATION FOR THREE 
DIMENSIONAL SUBSONIC FLOW BY A 
FIELD INTEGRAL METHOD 


& Thesis Submitted 
In Partial fulfilment of the Requirements 
for the Degree of | 


‘MASTER OF TEGHNOLOGY 


oe We “ ‘ as . s ‘ Z 
FF se gta ee koe ¥ ists , ex / 
‘ RO ek 
Se tue Sn 333 a 
* # 4 § % a i a. 2 
ae fe HS =, ‘ 4q : # . 


by 
KUMAR ALOK 


: | — to sane ) 
= "DEPARTMENT OF AERONAUTICAL ENGINEERING 


INDIAN INSTITUTE. OF TECHNOLOGY, KANPUR, . 
JANUARY, 1989 e 





vit, KAR PUR 


AE ding ALO- SOL 





oo 
ee 
my 


+ ten 


% . oe a 
anon. 





CERTIFICATE - REE reir pcs 


Certified that the thesis titled, "SOLUTION OF 
FULL POTENTIAL BQUATION FOR THE THREE DIMENSIONAL SUBSONIC 
FLOW BY A FIELD INTEGRAL METHOD" is a bonafide work carried 
out by Mr. KUMAR ALOK under my supervision. This work has 


not been submitted for a degree elsewhere. 





December 29, 1988. (Dr. N.L. ARORA) 


Professor ; 
Department of Aeronautical Enginee 
Indian Institute of Technology 
Kanpur 


174) 
ACKNOWLEDGEMENTS 


I wish to express my deep sense of gratitude to 
my thesis supervisor, Professor N.L. Arora, for his able 


guidance and constant encouragement at every stage of this 


work. 


ae would also like to sincerely thank Mr.J.P.Agrawal 
and Mr. Jyoti Swarup for their cooperation and help. I am 


thankful to Mr. U.S. Mishra for his excellent typing. 


December 29, 1988. 7 -~KUMAR ALOK 


CHAPTER 1 


CHAPTER: 2 


CHAPTER 3 


CHAPTER 4 
APPENDIX A.1 


A.2 
APPENDIX B 
APPENDIX C 
APPENDIX D 


APPENDIX E 
REFERENCES 


COMPUTER PROGRAM 


CONTENTS 


LIST OF FIGURES 
List OF SYMBOLS 
ABSTRACT 


INTRODUCTION 


MATHEMATICAL FORMULATION 


NOM DM NHN 
¢ 
mf G2 Ne 


The governing equation 
Integral equation formulation 
Internal Drichilet formulation 
Numerical discretization 


COMPUTATION SCHEME 


Sed 


Surface grid generation and panel 
geometry definition 

Field grid generation and volume 
elements definition 

Influence coefficient evaluation, 
assembly of coefficient matrices and 
iteration procedure 

Pressure coefficient evaluation 
Evaluation of tangential perturbation 
velocity component 

Central difference scheme 

Forward difference scheme 

Backward difference scheme 
Derivatives at the field points 


RESULTS AND DISCUSSION 


Evaluation of source influence 
coefficient 
Evaluation of doublet influence 
coefficient 


Evaluation of wake influence 
coefficient 


Evaluation of volume source 
influence coefficient 


Evaluation of the surface integrals 
over the singular surface go 


Introduction to Computer Program 


SAMPLE. INPUT AND OUTPUT 


iv) 


Vv) 
wit) 
xii) 


Lt 
3 
a7, 


20 
20 
23 


25 


44 


48 


ao 


he, 
64 


Number 


LIST OF FIGURES 
Title 


The definition of the three dimensional 
body surface and the region D 
Nomenclature for internal Drichilet 
formulation 

Definition of I and J lines in the surface 
grid 

Schematic of a surface grid 

A panel with its cormer points 
Quadrilateral y~ of the average plane 

3D grid in the field 

A volume element with its corner points 
Definition of geometric quantities for 
differentiation senate 

Points on the grid where forward/backward 


differencing is used. 


Definition of geometric quantities of a 


quadrilateral element 

Geometry of a general horse shoe vortex 
emerging from a wing Eraliing ecge 
Geometry of ote panels for a swept wing 
Exclusion of the singular point P from the 
body surface 

Cosine spacing 


Half cosine spacing 


v) 


66 


67 


67 


67 


67 


68 


68 


69 


70 


70 


Jen 


71 


72 


Te 


wee 


vi) 


Page 
Number Title 
. 73 
4.1 Comparison of chordwise C, distribution 
over a symmetric RAE wing 
4.2 Comparison of chordwise Ce ore 74 
over a lifting RAE wing 
4.3 Comparison of chordwise os Gistribution 75 
over a lifting M6 wing 
4.4 Comparison of chordwise Cy distribution 76 
over a symmetric rectangular wing 
4.5 Chordwise pressure distribution over 4 V7 
lifting rectangular wing 
78 


4.6 Chordwise distribution over a lifting 


rectangular wing. 


mi 


vii) 


LIST OF SYMBOLS 


free stream speed of sound 


Signed perpendicular distance of the 
projection of the point in the plane of 
quadrilateral from the extension of the 


side (Figure A.1) 


Ari? Ado? Ara 


Ao Anos 
Anne Anse 
oe 

32°433 


op 


A 
Aaj 


elements of transformation matrix 


any vector field as given by equation 2.8 


th 


potential induced by j volume source 


h control 


element of unit strength at i 
point 

wing span 

potential induced by gen source panel of 
unit strength at on control point 

the cosine of the angle between aa | 
Side of an element with respect to X axis 
defined by eqn. C.20a 


coefficient of pressure at - control 


point 


root chord 


h 


potential induced by ae doublet panel of 


unit strength at ioe control point 
potential induced by yuh wake panel of 
unit strength at oe control point 


inplane distances measured in differentiation 


o* 


n> D> 


Z 


viii) 


region exterior to the body surface 


Pigs. O24 5 


region enclosed by the body surface 


PAGS 9) 
unit vectors along t and y directions 


(Fig. 3.8% 


Function defining the body surface 
(3° 4 22172 


strength of fF Soins SsOurce element 


corner points of a volume Clement in 

global coordinate system 

the arc length measured along the side from 
the intersection of perpendicular with the 


extension of side i; the positive direction 


being from i to i+1 corner points (Fig. A.1) 


free stream Mach number 


unit normal vector 


unit normal vector grawn On the inside 
surface of the body 


number of volume elements 


Number of elements on the body surface in 


the flow 
Number of wake elements 


general field point where potential is 


evaluated 


average point of a surface element in 


global] coordinates 


Pie Poe Pe 


Va V ve 


x oO 
KX pV -Z 


ix) 


centroid of a volume element 

corner points of a surface element in 

global coordinate system 

integration point on the body surface 

i.e. the point where the singularity is 
placed 

distance from a PSG oink P to the 
integration point Q 

distance between the control point FP and the 
average point of the sukface panel or the 
centroid of the volume element whose influence 
is to be obtained at point P 

surface which encloses the region D 

surface of the body 

surface of the wake 


, side 


the sine of the angle between i* 
of an element with respect to X axis 
maximum diagonal of an element 
vectors in the plane of an element 
vector normal to the plane of an element 
normal and tangential velocities at 

7oe senerey point 

free stream velocity 

velocity components in global coordinate 


system at yuh 


control point 
assembled doublet and wake coefficient matrix 


Cartesian coordinates of a point in space 


XeVe 


TE 


R 


us 


vrr> 
3> 4 & 
® 
w> 


(no 


(7 Sy"s = 
Nn ah 2 


dimensionless Cartesian coordinates of a point 
in space 

X coordinate of the trailing edge in global 
coordinate system 

angle of attack 

ratio of specific heats 

doublet distribution on wake 

source distribution on surface panels 

polar angle | 

total potential at any point in the flow field D 
total potential on the inner surface of the body 
free stream velocity potential 

perturbation potential 

trailing edge sweep angle 

local sweep angle 


density of fluid 


o 3 o 
- ax 2 Wy j +3 k 
3 9 2 3 ; 
5 2 ee 2% 
x z 
surface far from Sp 


unit vectors .constituting a right handed 
orthogonal coordinate system attached with surface 


element 
components of &, nN, % respectively in global 


coordinate system 


chordwise and spanwise directions along 


which differentiation has been done 


first and second derivatives of @ with respect to 


X,Y,zZ crespectively 


xi) 


Pye Pye P, first derivatives of p with respect to X,YezZ 
respectively 

Superscript 

* quantities in the dimensional form 

Subscript 

CO free stream conditions 

CP surface control points 


gP field control points 


xii) 


ABSTRACT 


A three dimensional integral equation formulation 
for the solution of the full potential equation for subcritical 
transonic flow is presented. A general low order panel method 
with internal drichilet formulation is used to solve the 
incompressible flow problem governed by the Laplace equation. 
Strength of the piecewise constant surface source singularities 
is determined by imposing the external Neumann boundary condi- 
tion while doublet singularities distributed piecewise constant 
on the body surface is unknown of the problem. The body and 
wake surfaces are discretized into No surface and NU wake panels 
to obtain Ny, linear simultaneous algebraic equations in N, 
unknown doublet singularities. Gauss elimination technique 


is used to solve the system of simultaneous equations. 


The volume integral term of the Green's third 
identity is used to extend the incompressible flow solution 
to the subcritical transonic flow governed by the Pelisebn 
equation. A field grid in conjunction with the surface panels 
is used to solve the nonlinear equation. To do this, a finite 
parallelopiped shaped volume around the body is taken which 
is discretized into Ng volume elements where piecewise constant 
volume source singularities are distributed. The strength 
G of the volume sources is a function of total velocities bye 


a e.. and density derivatives p,, Pye p. etc. which is also 


unknown of the problem. An iterative scheme is employed 

to solve the full potential nonlinear ssGaeion. The conver- 
ged value of the velocity potential § is differentiated 
using second order accurate central difference scheme to 
Obtain the velocity components at all the surface control 
points. The isentropic flow relation is used to evaluate 


the coefficient of pressure C_ at these control points. The 


Pp 
results for chordwise pressure distribution over wings for 


both incompressible and subcritical flows are computed 


13,14,15 


and compared with the results of references The 


results show good agreement. 


CHAPTER 1 


INTRODUCTION 


The ultimate goal of computational fluid dynamics (CFD) is. 
to provide the aircraft designers with a method capable of 
solving the Navier Stokes equation over realistic aircraft 
configurations. Despite rapid advances in solution tech- 
niques and computing power available, this goal remains 

a somewhat distant target. A more tractable and thus 
popular approach is to consider solutions of the inviscid 
form of the flow equations and couple these with a boundary 


layer method. 


A great deal of progress has been made in solving 


linear and nonlinear inviscid flow problems using finite 
difference and finite element methods (field methods), 

but these methods have not yet proved to be easily adaptable 
to complex three dimensional surfaces. The field methods 
are capable of solving the full potential or Euler equations 
but with a more limited geometric flexibility at present. 
This limited flexibility of field methods arises from the 
requirement to use a computational space or grid which 
eve: qs region sxtechal to the eee but which 


must also conform to the configuration surfaces. However, 


it is a very complex problem for a general three-dimensional 


geometry. 


The best developed and most widely used methods 
for computing flow about arbitrary geometries are those 
which represent the configuration boundary with some form 
of source/doublet-vorticity distributions, generally called 
surface panel methods. Panel methods do not have the 
geometric problem as in the case of field methods, because 
surface definition only is required. Historically panel 
methods have been limited to solving problems governed by 
the linear Prandtl-Glauret equation. Recent studies have 
shown that supeace panes used in conjunction with field 
grids can also solve nonlinear problems. Besides their easy 
adaptability to complex aircraft configurations, these 


methods provide nearly exact results. 


The first work on the flow field around three 
dimensional bodies using panel methods is the one by Hess 
and smith. Their method uses constant strength source 
elements to solve the problems of steady subsonic flow 
around non-lifting bodies. The strength of the source 
@lement is obtained by imposing the tangency boundary 
condition. This method has been extended to lifting bodies 
by Rubbert and Saaris*. They have used doublet singularities 
in addition to source singularities. Panel methods have 


been successfully extended to solve unsteady subsonic and 


a 
supersonic flows also, as shown by Morino . 


There are various schemes of panel methods which 


can mainly differ in the following ways: 


1. 


3. 


The choice of surfaces on which sources and/or 
doublets are located i.e. those schemes may 
consist of the wetted surface onlyor fictitious 
internal surface only or both of these. 

The choice of unknown functions on these surfaces, 
such as source density or doublet density or 
vorticity or a mixture of these. 

The form of the boundary condition such as 
internal or external Neumann and internal or external | 
Drichilet conditions. In a Neumann boundary 
condition we specify 36/3n (i.e. normal deriva- 
tive of velocity potential), while in a Drichilet 
condition we specify the velocity potential $ 
over the surface of the body. <A condition is said 
to be internal or external depending upon the 

body surface (i.e. inner or outer) where it has 


been satisfied. 


Another classification is based on the type of poly- 


nomial which approximates the singularities on the panels i.e. 


low order and higher order panel methods. In low order panel 


methods piecewise constant surface singularity distributions 


are used. 


These methods offer great flexibility for appli- 


cation to general problems.Since continuity of singularity 


distributions is not enforced from panel to panel, the 
assembly of panels representing a complicated surface is 
fairly straightforward. This method has been discussed 
by Morino® and Maskew?. The low computing cost makes it 
practical to apply the method to nonlinear problems requiring 
iterative solutions e.g. wake relaxation for high lift and/ 
or multiple component problems, viscous/inviscid calculations 
with coupled boundary layer codes including the case with 
extensive separation: © and also eine stepping calculations 
for three dimensional unsteady problems*’, Since in low 
order panel weeKeas singularities are not continuous at the 
edges there is some leakage of the flow, to reduce this error 
higher order panel methods have been employed by Johnson*®, 
In higher order panel methods we have more continuous repre- 
sentation of the surface singularity distributions. With 
this type of representation these methods should allow a 
reduction in panel density for a given solution accuracy, 
and hence should lead to lower computing costs. No such 
benefits have appeared so far for the general ieee dimen= 
sional case. But in the sensitive flow regions ewes near 
the wing leading and trailing edges and wing flap attachment 
etc.) higher order methods may be more reliable. 

As mentioned earlier, panel methods in conjun- 
ction with field grids can solve nonlinear problems +? >¢ 6+ 7 
also. Since panel methods are usually considered to be 


superposition methods, their extension to nonlinear problems 


may seem surprising. The theoretical basis for this 
extension is shown in Chapter 2 using Green's third iden- 
tity?. BaD Jonnson> and Larry L. Erickson® have solved 

the nonlinear equation using a related but more powerful 
method. This method is iterative and each iteration 

requires solving two poisson equations which are obtained 
from the optimal control approach (i.e. variational calculus). 
The poisson equations are solved using surface panels and 

a grid of field points. The panels are used to enforce 

the boundary conditions, and the field grid is used to 
evaluate the volume integral by fast Fourier transforms 
(FFTs). The method described by sinclair’ finds the velocity 
induced by the surface and field singularities at all the 
surface and field control points to solve the full potential 
nonlinear equation. He Pivse.goiyes the Prandtl Glauret 
equation to obtain the velocity components at all the surface 
control points-by approximating the surface geometry as a 
set. Ge Planar quadrilaterals and placing source and vorticity 
singularities at these quadrilaterals. A region of space 
including the configuration and its immediate neighbourhood 
is taken to extend the Prandtl Glauret solution to transonic 
flows. This space is subdivided into cuboids where field 
sources are distributed. An iterative scheme is employed to 


solve the full potential Poisson equation. 


The present report is aimed at solving the full 


| potential equation for subcritical transonic flow using a 


low order panel method. An internal Drichilet formulation?” 


is used to solve the incompressible flow governed by Laplace 
equation. The unknown of the problem are doublet singulari- 
ties piecewise constant at the panels. Source singularities 
are known from the external Neumann boundary condition. The 
whole surface is divided into Np panels where piecewise 
constant source and doublet singularities are distributed. 

On discretization we get N, linear simultaneous algebraic 
equations in terms of unknown doublet values at the XN, panels. 
Gauss elimination technique is used to solve the simultaneous 
equations. A second order accurate central difference scheme 


is used to find the velocities at the panel control points. 


To solve the full potential nonlinear equation, a 
finite volume around the body is taken which is divided into 
Ne cubic elements where piecewise constant volume source 
Singularities are distributed. The nonlinear flow is governed 
by Poisson equation, i.e. oo = G, where G is a function of 


2 Pie Gs 2608s etc. The incompressible flow 


x * Pxx Pye Pyy 


solution technique has been extended to compressible nonlinear | 


ene (Chapter 2). 


problem as shown by Green's third identity 
The incompressible flow solution is used to find ¢§ at the 
centroid of the field elements. Second order accurate central 
difference scheme is used to evaluate various derivatives of 

g@ and hence source strength G is evaluated at all the field 


elements. The influence of these volume sources is evaluated 


on the surface control points and we solve the modified 


system of linear equations in Ny unknown doublet singula- 
rities. Again the influence of surface singularities is 
evaluated at the field control points. This iteration 
procedure continues until we get a converged value of $4 

at all the surface and field points. This method converges 


very fast and gives quite accurate results. 


The mathematical formulation of the problem and 
numerical discretization is presented in Chapter 2. The 
computation scheme used to solve the problem is explained 
in the Chapter 3. In Chapter 4, the results computed 
over three types of wing configurations are presented and 


13,14,15 — The closed 


compared with the results of references 
form expressions for surface source and doublet singularities, 
wake surfaces and field source influence coefficients are 

presented in Appendices A,B and Cc. Evaluation of the surface 
integrals over the Singular surface g is given in Appendix D. 


Finally an introduction to the Computer Program is presented 


in the Appendix E,. 


CHAPTER 2 


MATHEMATICAL FORMULATION 


2-1 The Governing Equation 


If a three dimensional body is immersed in a steady 


inviscid, irrotational flow, there exists a velocity potential 


* 
@ which satisfies the following equations: 


Conservation of mass 


Conservation of energy 


2 2 2 2 
x? =_j ..* _ * Y -1 ao x 
a + a V = a ae ( ¢$ + 2 + 


Isentropic gas equation 


* * a * * ¥ ~1 
p/P, = fa fay) 


* * 
where 6 is the total velocity potential, a the 


(2.1a) 
2 
* 
o ) 
z 
(2:1) 
(2.1c) 
speed 


* 
of sound, p the density,y the ratio of specific heats, 


| x 
and subscript oo for the free stream conditions. We use Vi60? 


| x | 
P and characteristic length c (of the body) as reference 


parameters to obtain the dimensionless form of the above 


equations: 


(2<id) 
* * * *x 
P= PSP MO = VG AG 
Combining Eqs. (2.1) we get the following dimensionless 
equation: 
2 
b = _ * 
xx +? Oy oP G or gy $ G (2e2aj 
where 
1 | 
= - = (Pp a 
G 6 ey o. es ee ee (2.2b) 
and 4 
= Y -1 2 2 _g2 _ 32 ¥ =1 
SS a ee ML Mae o., a) (2.2c) 


If the free stream Mach Number is such that there is no 
formation of shock wave anywhere in the flow field then $ 
is continuously differentiable and has continuous partial 


derivatives of the second order in the entire flow field D. 


¢g and its derivatives must satisfy the following 
conditions: 
(i) ¢ must satisfy the external Neumann boundary 
condition at the surface of the body, i.e. normal 
derivative of the total velocity potential $ 
should be zero at the body surface. Let the 
body surface be defined oy dhe equation 


F(x,yez) = O | (252) 


1-9 


The above condition may be stated as 


2e. ra ae 98 |  G (0.4) 
F=0 


where 1 is the unit outward normal vector at the body surface 


and is written as 


hRo= (9 F/1 VF! )gLg | C2250 


(4i) & approaches free stream value of velocity potential 


at infinity 
$—e $ at infinity | (2.6) 


Total potential ¢$ can be written as a sum of free 
stream potential 6, and perturbation potential @ (due 


to the presence of body), i.e. 


6 = § + @ (2.7a) 


CO 
where 
v2. os = Yoo (2.7b) 

and 

seg 

Vig = 0 (216) 

ae gg = 9° (8 +9) =v°9 = G (2.7d) 
where 


G is given by Eq. (2.2b). 


a2 


2e2 Integral Equation Formul ation 


Using Green's identities we can express @¢ in 
terms of singularities distributed on the surface of the 
body and in the entire flow field. First we'll establish 
the various Green's identities°’?? and then use them to 


solve Eq. (2.2a). 


For a vector field A, the Gauss theorem gives 


ds = - L[{[f[IVA AV (223) 
V 


where ‘XN is the unit normal taken +ve into the volume D. 5 
consists of all the surfaces that bound the region D, see 


Fig. 2.1. 


Let $d and ¢$!' be two continuous functions with 
finite, single valued, first and second derivatives continuous 


throughout the volume D. 


“We define A = 696! (2.9) 
we <- t 
AR = afg @'. RP = & B- (2.10) 
yA = 9. ($9%') = 9S" 455.58! a5 


substituting Eqs. (2.10) and (2.11)in (2.8) we get 
3g ' 7 | 2 
fi on ds = - fff [vee $! +O¢9 6! } av 
S D 
| (2.12) 


Equation (2.12) is Green's first identity. 


LZ 


Interchanging $ and » ' in Eq. (2.12) we get 


~ x 2 
i gt oo ds = - siipg'.gete' gs Jav 
D 
(2213) 
Subtracting Eqs (2.13) from Eqy (2.12) 
ff og : oe oe ' 2 | 
S en ds 7 aa an ds ieee Se ee 
~ IT § 9% sav (2494) 


Equation (2.14) is Green's second identity. 


Let 6 be the velocity potential which satisfies 


the Poisson eqn. es = G and 


Pg * Sccy, ject ee ee (2215) 
Y (xp-%, a (Yo-¥p)* + (Zp-Z : 

where r is the distance between a field point P(x), Yp zo) 
in the region D and Q(x, For Zo) the variable of integra- 
tion. p | satisfies laplace eqn. ee = 0 everywhere 
except in the near vicinity of point P where §' has an 
impulsive behaviour. So we exclude this point from D. 
This we do by centering a small spherical surface 4g 
around the point P. Now in the resulting volume DY 3 = 


is harmonic (i.e. it satisfies Laplace equation) everywhere 


and we have 


13 


fi O 1 1: op a 1 
6$— (=) ds - Jf = — ds + ff5— (=) do 
S-o en = S-o on Oo on 2x 
~ sf 2 a esl 4 v 8 av 
Oo D* 


(2.16) 


If we let go become a very small sphere then (See Appendix D) 


Q 1 | 
lim sf/6 SS (2) ds = -4n @ (2.17a) 
o 2 0 
and 
= : 
lim ff = 5m 4s = 0 (2.17b) 
o-0 


Substituting Eq. (2.17) in Eq. (2.16) we get 


fr Be . 


This ie Green's third identity-The first term on the right 

is the potential of a doublet distribution on S of strength 
(@/4n), the second is potential of a source distribution 

of strength ( - $50 8/ On) while the third is the potential 
of a wolume source distribution of strength - v8 / 4n ° 


2-3 Internal Drichilet Formulation 
Let $ and o, be the velocity potential in 
the region D and D; respectively which are separated by 


the common boundary S,,(Fig.2.2). D is the external flow region 


14 


enclosed by S, (the body surface), S,, the wake of zero 
thickness across which + is continuous but $ is dis- 
continuous, and © the outer boundary extending upto infinity. 
D; is the internal flow region enclosed by Se See Fig. 2.2. 


Let P (x,y,z) be a point in the region D,using Green's third 


identity (Eq. 2.18) for region D we get potential at P as 


a Ls ey 1 re a Gd 
P . On oF ot Ae (ds 
47 Sats on 4 S.+S. +7 an 
BW 
= fif i 


Since fluid is at rest at infinity (i.e.¢ is bounded) 


integral over yr yields Pes 


=) os = == ae We) 
P oO 47 Sptsy on ‘Ir 47 SptS,, gn r) = 
ae Lee 
ae [Lf = v< ¢ av 
D 
(2.20) 
Since P lies outside D, so = is harmonic everywhere in 
region D; (i.e. it satisfies Laplace eq. is * = Q) 


everywhere in D;). Using Green's third identity for Ds 


we get 
1 rs 3 1 1 a 1 
0 = = @# — (=) ds - = 44 — (=) ds 
47 S 5 1 on; = 41 S an, © 
an Dd; . VG, av, (2277) 


i 


at any point on S5,n and - are in the opposite direction, 
3 a 
so an; = = an 2392) 


adding (2.20) and (2.21) and making use of (2.22) we get 


a 1 3 
$ = 1. ff (6. 6 sa, = 
: a r= $5 ( ) a (=) ds 
ag 
- = —_—> - — S 
a 
1 FSS 250 2 1 a. Peo fh 
~ 4a D+D, (ye ~ ¥ oH avo + Sy oan (>) dS 4 
t. hf 22 pe 
iz 8! gn (=) a5, (2523) 
| og ; 
last term in Eq. (2.23) is zero because a is continuous 
across the wake. Since $ is discontinuous across the wake 
and n,=- ny, we get 
cof 8. gle oe _ so ye 
26. $—5 () asy = OG - FL) Se OC) as, 
Sw u 
(2.24) 


substituting Eq. (2.24) in Eq. (2.23) 


. 36 | 
= 2 SF ge ol ~ Loge Of . ay as 
e, = c.. & Be S (g- 2,) aa (=) ds es a os =). 05 
5 
2 oe ee SLI 2,264 1 i 
oa 3, (o,-?,) on, (=) ds + D+D, C9 $-9 ¢,) = av 


(2425) 
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If P lies on S, then we'll have a factor of 2n instead of ce 


B 
Proof of this is given in Appendix D. 

Equation (2.25) will produce same potential $6 p for 
various combinations of $ and ? ye 060 For uniqueness of the 
solution we must specify either 6 or $ 3 or their normal deriva- 
tives. This we do by specifying $ 4 =? in D, to ensure 


oO 
that there is no perturbation inside the body. 


We have 
SS 2. oF (2. 26a) 
$- 6, = db. $F 8 (2.26b) 
ag 8 oto 8g 

9 an gn an * “On. = On (2.26c) 

2 2 _ 2 2 2 x; wie 

yv@ ~V¥ %, = WV &, + VO VY os = 9 
(2. 26d) 


Using (2.26) in (2.25) we get 





2 a a9 1 
Ep 9p = ff ea Gl as- 44 se os as 
Sp S , 
Bo 
+ Jf @ -9,) ; (4) ds : ffflog 7 P? av (2 27) 
W my 4 3 my 7 D+D, 
where En = Ay, if P lies in D but not on Sp 


= 2m if P lies on S, 


if P lies in Dis ?, = 0 


(i) 


OL, 


(ii) 


7 


Value of a2 is known from the boundary condition 


on 
(2.4) 
3g 
ao 
a /, 7° @ OER t+ Hs 0 
F=0 
ce = - of a a, ni 
on VP x ° Gp * 
98 LY a (2.08) 
an Cc ° | ° 


Strength of the vortex sheet emanating from the trailing 


edge is taken to be the value of @ at the centroid af 


the adjacent panel i.e. 9, -P, =APoe > This is 


how we satisfy the Kutta condition at the trailing 


edge. 


wake 


There is no variation in the strength of the 


in the streamwise direction but has a variation 


in the spanwise direction. 


2e4 Numerical Discretization 


The integral Eq. (2.27) is discretized in the 


following manner: 


1. 


Surface Sp is discretized into flat quadrilateral 

e — { - > | 2g 
panels £, 7; k = 1, Np. Values of g and Vege? ( 5m) 
in each element are represented by those at its 


average point and assumed constant in each element. 
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2% Surface Si is discretized into semi-infinite 
strips D1, w=i1, Nj. Value ofap= 9, - ?, is 
assumed constant in each strip element. Values of 
Py and P45 are approximated by those of o at the 
panels adjacent to the wake strip on Sp: 

a | A finite parallelopiped shaped volume is taken 
around the body for volume D+D,. This parallelopiped 
is subdivided into small cubic elements V,; 

k = 1 to Noe Values of » and G in each element 

are represented by those at its centroid and 


assumed constant in each element. | 


If we choose average point of the os panel on 


. the body as our control point then Eq. (2.27) can be written 


as 
a y | 
2m 9; = iz , Coy + E GO, Bi, + oa Mw es 
j#¥i j=1 
N 
g 
+ a Gy Bog 
gg 
i=1, Np (2.30) 


similarly if we choose centroid of the had volume element 


as the control point in the field then we have 


Ae - a | Ng 
4no, = 2% 9, C., + 4 GF BL. + : C + 3S Cx Be 
a i j=1 J 13 - j=l j 1) j=1 b Ma W434 cn J ij 
foe 4.413. (2.31) 
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where Cc. GC and Aa; are the potentials induced 


4 ax cas 

Se aa Wag 

due to the unit strength of as doublet panel, source 
Panel, wake panel and volume source respectively at 


the -? control point 


0 1 
and c.. = f —— (——) dz. (2.32a) 
+J t, any. Fay 
J 
B 4 i e (ntmemcanem e32 
ij ae cs ay, (2.32b) 
j 
oO : 
om zs = yffo os (=) ax. /(2.32c) 
ij Ws u, ij 4 
_ i. | 
As ee FIL ao. dv (2.32a) 
V 1j 
j 
"Ws = (P., - Fy), (2.33a) 
NN A | 
— i. 7 
Gj ~ P Pa be, BP oP 2): (2.33) 
ify-l 
= Yes . 2 aay 2 ae 2 oe aC ‘ an 
P, [i+ pM (2 oy  8y - 22)7; (2.33d) 


Closed form value of the coefficients Bij and Cie ome and 
ij 


Aj; are given in Appendices A,B and C respectively. 


CHAPTER 3 


COMPUTATION SCHEME 


3.1 Surface Grid Generation and Panel Geometry Definition 


To solve the flow problem, the body surface is 
discretized into a number of flat quadrilateral panels. At 
the edge, the panels have continuous position but discon- 
tinuous slope and curvature. The panels are so closely 
spaced that they do not deviate significantly from their 


average planes. 


To generate the nanete, the wing surface is 

divided into a grid of corner points Pit, as I= 1,M and 

J = 1,N where the components of the position vector P are 

in global coordinates (x,y,z). The index I stands for 
lines of constant x/e (i.e. x) and J for lines of constant 
2Y¥/b (i.e. YY). Index I increases from leading edge to 
trailing edge along the upper surface and from trailing edge 
to leading edge along the lower surface. Index J increases 
from midspan to the tip of the wing. This is shown in 

Fig. (3.1) . The M lines and N lines in the chordwise 

and spanwise directions respectively can be spaced either 
uniformly or by some other law to enable us to have more no, 
of panels near the wing leading and trailing edges and also 
at the wing tip. Since relatively large perturbations occur 


near flow regions mentioned above, the latter choice is 


preferred and we use cosine and half cosine spacings in 
the chordwise and spanwise directions respectively. 
Schematic of the surface grid is shown in Fig. (3.2). 

A panel is defined using four adjacent corner points 
P(I,J), P(I + ino, (PCr, J+1), P(I, J+1) and are 
relabelled as Box 36 P,and ey respectively. Panels are 
numbered such that for a constant J index, panel number 


N increases with Index I. 


iee. N = (J-1) * (N=1) + Tp (J = 1, N-13 I = 1, M-1) 


A panel with its cormer points is shown I Figs 3236 


It is convenient to evaluate the source and doublet 
influence coefficients in the panel coordinate system. Hence, 
we define a local orthogonal coordinate system attached with 
the panel. Let T, and T., be two vectors lying in the plane 


of the panel where 
a = Pp ee Px. a P (3.1a) 


and 


34 
if 
Os 
4: 
ig 
| 
ty + 
{ 
A) 


7 3 4 (3.1b) 


Outward normal to the plane T is defined as 


3 


tT. = T, x T, (3.1c) 


cap _ = 24 . | | 
oe oS eo, (3.14) 


dd 


vectors t T, and c are mutually perpendicular and their 


a * 
unit vectors es 1 e z respectively define the right handed 


| orthogonal system attached to the plane of the panel 


22 











a 
A z 
— 3.2a 
és iT | ( ) 
| tT. 
§ 
ITs | 
oh 
> = = : | (3426) 
[TI 
The average point P, of the panel is defined as 
> as 1 = > at a | 
Boo Gly, Py Se Ry ee Be (3.3) 


The plane passing through P with normal * is 


defined as the average plane of the points Py er P, 

and B,. This plane is average in the sense that it passes 
through the midpoints of the line segments joining P, 

and Po, B, and Pas P, and Py and Py and P, as shown in 
Fig. (3.4). The quadrilateral 5 thus formed is called 
average plane. The local coordinate system associated 


with 5 consists of three orthogonal vectors ey De % 


with 
origin at the average point Po: The average point P, is 
also taken as the control point of the panel to satisfy 


Eq. (2.30). 


With this choice of local coordinates any point 
P with components (x,y,z) in global coordinates can be 
transformed to the point O(& n,2Z ) in the local coordinate 


system using the transformation 


0 = A. (pie Py . (3.4) 
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i.e. 
"ot, “a7 AL3 | os 


il 


21 22 23 | 


31. 432 a3 | eat < 


where A is a transformation matrix whose component are 
the direction cosines of oe ae with respect to global 


Nn 
coordinates (x, yrz)s ; a and z% are defined by Eq. (3.2), 


A A A 
Sy Sy — 
- A ~ “A 
fA] = lx Ny Ne ! (3.6) 
A nw “A | 
oi oy S | 
ad 


3.2 Field Grid Generation and Volume Elements Definition 

As we have seen earlier in Chapter 2, there are 
volume source singularities distributed in the flow field to 
represent the compressibility effect. We construct a finite 
parallelopiped shaped volume around the wing. This volume 
extends upto half M.A.C. (Mean Aerodynamic Chord) ahead 
of the leading edge at the root and half M.A.C. behind the 
trailing edge of the wing tip along x axis, one M.A.C. above 
the XY plane i.e. along the z axis and 1.5 times half span 
along y axis. This volume is subdivided into a number 


of cubic elements. 


24 


First of all we generate a 3D grid inside the 
parallelopiped of cormer points G(I,J,K); I=1, Nx; J=1,NY 
and k=1, Nz; where the components of position vector G are 
in the global coordinates (x,y,z). The indices I,J and K 
are for the lines of constant X,Y and Z and increase along 


the respective directions. 


Since flow is more perturbed near the wing leading 
and trailing edges, wing tip and near the wing surface in 
general, we prefer to have more number of cubic elements of 
variable sizes in these regions. See Figs. 3.5a and 3.5b. 
Hence, along the x axis we use half cosine spacing in the 
region ahead of the leading edge and behind the trailing 
edge and cosine spacing on the chord, along y we have half 
cosine spacing while along z the Nz lines are spaced in 
the ratio 1:23:3....: (Nz-1). Alternatively, we can have. 
cuboids of uniform size distributed all over the flow field, 
but this choice is not capable of capturing the perturbations 


properly unless we take avery large number of cuboids. 


Using the eight adjacént corner points G we define 
acubic element (see Fig. 3.6). The corner points of a 


cuboid are numbered as: 


CG.ucn) 2 G, > G(I,d,Kti) = C. 
G(It1,Jek) = Gy ; G(I41,J, KH) = 6 ied 
G (I41, J+1, K) = G, >; G(I+1,741,K+1) = e 
G(I,J+1,K) = G, ; Cie = Ge 
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To number these elements, for a fixed Z plane first 
we move along y then along x direction. After numbering this 


2 plane we number the next Z plane in the similar manner. 


in other words 


N 


Hi 


(K-1) * (Nx-1)* (Ny-1) + (I-1)* (Ny-1) + JO; 
K 


i 


1, Nz-1l; J = 1, Ny-1;7 I = 1, Nx-1 (3¢6)) 


The volume coefficients are evaluated in the global coordi- 
nate system, so we don't need to define any local system 


attached with the cuboid. 


The centroid PB. of the cuboid is defined as 


= 
Po = 


Olr 


[ G, + G. + Gy + Ga + Ge + Go + Go + Gp j (3.9) 


The centroid Pa is used as the control point of the cuboid to 


Satisfy eqn. 2.31. 


3.3 Influence Coefficients Evaluation and Assembly of 
Matrices 


After defining the geometry of the surface panels 
and volume elements we evaluate the foi towing influence 
coefficient matrices using the expressions given in 
Appendices A,B and C. 


(1) C. : influence of doublet singularities on the surface 
P control points 


(ii) C__s influence of doublet singularities on the field 
JP control points | 


(iii) B__: influence of source singularities on the surface 
cP control points : 
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(iv) Bas : influence of source singularities on the 
field control points 

(v) C+ influence of wake singularities on the 
gp field control points 

(vi) Cc, + influence of wake singularities on the 
cp surface control points | 

(viilj Aap : influence of volume source singularities 

| on the surface control points 
(viii) Pes : influence of volume source singularities 


on the field control points. 


if the flow is symmetric about the Y=0 plane, we 
will have same perturbations at the points which are mirror 
images of each other about the Y=0 plane. Hence, strengths 
of the corresponding surface or volume singularities at 
these points will be same; so we can add the corresponding 
“in€luence coefficients. Again we can add the doublet and 
wake coefficient matrices because the strength of the wake 
surface is assumed to be the same as that of the panel 
adjacent to it on the body surface, and if a panel is not in 
contact with the wake, the corresponding wake coefficient is 


set to be zero. Let the assembled wake doublet coefficient 


matrix be W. Equation (2.30) can be written as 
[Wan] {> op} = [Bop] eae ae [Aap] 1Gi 
Or | 
N N_x1 N x1 N_ x1 
"pew p" p Pp eeaa 
Wes! (Gag) [bo } + ff.>} 


similarly Eq. (2.31) can be written as 
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{ 7 | 1 
Rp 3 L Wop ee ae (o} +[ A] (6) 
or 
N_xl | 
a AS ae ea A 
Lo so = [Wend {? od + tbo, } iti 


(3542) 


Integral Eq. (2.27) has to be iteratively solved because ¢ 


and G are both are unknownse 


Iteration Procedure 
Ls First we make an initial guess for G (i.e. G=0) and 
solve Eq.(3.10) for £@ gph on the body surface using 
Gauss elimination technique. 
2s Knowing {9 oo we find Pap! at the field points 
using Eq. (3.11). 
3. Once {o gp} is known, we find G at all the field 
| elements using Eqs. €2«33¢) sand. (2.33d) 
4. This G value is used to solve Eq. C3610) 
36 Steps 2-4 are repeated until convergence occurs. If 
9 values in two successive iterations at all the 
control points on the surface and in the field differ 
by a value less than 107° then the @ values are 


supposed to have converged and the iteration procedure 


is stopped. 
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3.4 Pressure Coefficient Evaluation 

The values of perturbation potential @ at all 
the surface and field control points are thus determined. 
We evaluate the pressure coefficient Co at all the control 


points on the surface using the isentropic relation: 


2 
TI V a1 
2 ae? 2 ) vai 
Ee a See a eee IM ee ee) - 1} 
i Y Moo | : ie ns 
(3.12) 


oe 
where Vs is the net velocity at the qe8 control point. 


Ve. = V + i + v2 (3.14) 


2 2 2 a: 2 
Vi = VE + Vy =F VE. 3615) 
a 2 i. 
where Vy = O because of no penetration boundary condition. 


i 


3.5 Evaluation of Tangential Perturbation Velocity 
Component 


We choose four neighbouring panels surrounding | 
the panel on whose control point the tangential velocity 
and pressure coefficient are to be evaluated (See. Fig.3.7a). 
Using second order accurate central difference scheme along 


nN, and A, directions we evaluate a9/@ A, and ag/ or > 
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respectively. On the boundary of the wing (i.e. panels 
near the leading and trailing edges, root or tip), second 
order accurate forward/backward differentiation scheme is 
used. Such panels are shown in Fig.(3.8, F and F x 
represent those points on the body surface dea ne forward 
differencing is used to evaluate — and og/ ar. 
respectively. Similarly te and B k represent those 


points where backward differencing is used to evaluate the 


derivatives in NG and Ao directions respectively. 


O 
Finally = and oe can be resolved in two 


orthogonal directions in the foflowing manner”; 


a ° a2 | (3.16) 
1. 
De, ae es a a 6 VBA G.e 
ay = ak. 2° ©, . My ge) (3.17a) 
2 
or 
= = f a - Sing 99/3 / Cos Q (332-7) 
1 


where g is local sweep angle. 
“Net velocity at the control point in the t and y 


directions for unit onset velocity is 


3 
Vv, = Cos@ + =e (3.18 a) 
| a 29 | (3.18b) 
ee ay ; 
ee: 
a = (Cos39 + ao? + ie (3.19) 


dt oY 
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where @ is the angle that the onset flow makes with the 


panel in question. 


3.6 Central Difference Scheme 
With the help of taylor series expansion about 
point O we can approximate the value of at points 1 


and 2 (Fig. 3.7c) by 


2 2 
% ds 9 Po on 
Ps as ?, 7 9, ee + aot ax | + O(BA;) 
| - (3.908) 
ey at 37 
ss O° 2 Po 3 
~*% ot oe Fi ee 
(3.20b) 


multiply (3.20a) by d* and (3.20b) by (-a*) and add them 








2 
2 2 a 2 2 ° 54 
d5 91-4, Py = PY (dQ - Ay) - 4,4, (4, +4.) aK ; + O(A)) 
(3292) 
| = 20 2 
é 9 (r°-1) + 9, - rg 
Fo | = — Z : + 0 (AAS) (3.22) 
Ory Io a, (1+r) | 
where 


Here the truncation error is of the order of 0 (A *) which 
gives a very accurate approximation of the derivative 
when panels are closely spaced. Similarly using (3.22) 


Q gf er, can be evaluated. 


ae 


347 Forward Difference Scheme (Second Order Accurate) 


This scheme is used for panels near the leading 


edge and root. We evaluate $, and Po using Taylor series 


expansion about point 0 (See Fig. 3.74) 


i ap ag 
, =o +4 x f. 2S - o (ax? 
1 —— eOGx) “CHc2355) 
O 1 on, é 2 arnS 25 2 
; : : oi ( as 52 Re ands 
= + = eee ey GS ) 
2 O 2. Oka 1. 3 ane Lo 1 
(3.23b) 
, ies 
multiply  (3.23a) by a5 and (3.23b) by (-d})‘and add 
them to obtgin 
9, a2 -9,a* = 9. (a2 ~ a2) + (a, a2 - a, a2} — a! 
po a4 o > 1 1, De aha ie 
+ 0 (AA?) (3.24) 
or 
2 2 
| 9 6c -~ 9 =~ Pf (ro = 1) 
a 1 2 29 + O(AAZ) 
ony Oo a 2 
2 (r= 4) | 
(3.25) 
where re = d./d, - 


3.8 Backward Difference Scheme (Second Order Accurate) 
This scheme is used for panels near the trailing 


edge and the wing tip. 
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Again using Taylor series expansion about O(See Fig. 3.7e) 








z 2 
qd a 
a ee ae es eee Oe) 
1 O- 1 ONG 2 aX | 1 
1S 
(3. 26a) 
as oo 2 p ace 
P = ae - d. ah + d5/2 O gPA,| + a 
© 
(35265) 
multiply  (3.26a) by d$ and (3.26b) by (-d*) and add 
them to get | 
2 a See 2 2 2, 
Py 4g — Pe Fy = PQ bdy - Gy) + (dy ay - 4,45) sy 
: <s : “a, 20 
AnS 
+ O¢ 1) e727) 
2 2 
aC | , F = Go - F (oo a) | 
st Po ene eee eeieeneea (3.28) 
Lio d., (1 - x) 
where 


3.9 Derivatives at the Field Points 


x! 


We need to evaluate derivatives like 9@,, Pye Cos p 


Pi. PP. etc. at the field points in order to compute the 


Y Z 
volume source strength G. We evaluate these derivatives 
using central difference scheme as mentioned earlier in 


the Chapter. If the field points are adjacent to the 
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sides of the finite parallelopiped volume D then forward/ 


backward differencing is used. 


In the field perturbation velocities Py Pye Py 
_are along three orthogonal directions x,y,z. The net 


velocity along x,y,z direction is 


V, = Cosa + 9, | (3.29a) 
Wes = Py (3.29bD) 
V, = sina + 9, (3..29¢)) 
vs = ve + ve + ve (3. 29d) 


where q is the angle of attack and Vy the net velocity 


so.a° control point. Using eq. (3.12) Cp can be evaluated. 


CHAPTER 4 


RESULT AND DISCUSSION 





Three different types of test configurations are 
chosen for comparison of chordwise pressure distribution. 
The first configuration is a swept, tapered RAE wing which 
has a midchord sweep of 3075 aspect ratio 6.0, taper 
ratio 1/3 and NACA 0015 airfoil section. The second one 
is an ONERA M6 wing with leading edge sweep 20; aspect 
ratio 3.86 and taper ratio 0.583. The third and last 
configuration is a rectangular wing with aspect ratio 6.0 
and NACA 0012 airfoil section. All the three configurations 


are without camber,twist or dihedral. 


In all the computations mentioned below, surface 
Panels are spaced according to cosine law sas the X dire- 
ction and half cosine law along the y direction. The field 
elements £01 low half cosine spacing ahead of the leading 
edge. and behind the trailing edge and cosine spacing on the 
root chord along the X direction, half cosine spacing along 

the y direction and along z they are spaced in the 

ratio 1:2:3:.: (N-1) to have more number of elements near 


the regions where large perturbations occur.- 


Solutions from the method of Roberts’ * are selected 


as datum for all RAE wing cases. Only incompressible flow 
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equation has been solved for RAE wing cases. Solution from 


14 


the method of SALMOND is used as datum for M6 wing case 


while solution of se1it> over 2D airfoil is compared with 
the solution over a rectangular wing at its root, because 
the flow behaves essentially 2 dimensional near the root 
of alarge aspect ratio rectangular wing. For the last 
two test configurations we have solved the full potential 
equation for subcritical transonic flow. All the flow 


calculations have been made on Nexus 3000 Apollo system. 


tndcmpreseiote: nrow Solution 
For RAE wing configuration we have solved only the 

incompressible flow equation. Figures 4.1 and 4.2 show 
chordwise tessuve Bist eimucien On a RAE wing with NACA 0015 
airfoil section at 0° and 5° incidence respectively at 

i. = 0.625. The results are computed using 32x4 surface panels 
and compared with Roberts’ * datum solutions which use 42x13 
surface panels at n = 0.549. For both the cases results 


show very good agreement with the datum solutions. 


Full Potential Equation Solutions 


Full potential equation has been solved for ONERA M6 
wing and rectangular wing configurations. Figure 4.3 shows 
chordwise Cy, distribution on an M6 wing at free stream Mach 
number 0.5 and angle of attack 3°. This result is computed 
using 30x4 surface panels and 25x5x8 field elements and 


compared with Salmond-* datum solution at n = 0.659 which 
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uses 65x9x21 field grid. The solution converged in four 


iterations and is in good agreement with the datum 


solution. 


Figure 4.4 shows C., distribution at the root of 
a rectangular wing with as NACA 0012 airfoil section and 
aspect ratio 6.0 at 0° incidence and free stream Mach 
No. 0.72. 32x4 surface panels and 25x5x8 field element 
are used to compute the result which is compared with 
the solution of Sell over airfoil section. Figures 4.5 and 
4.6 show Cy distribution at the root of the same rectangular 


Oo O 
wing at Moo = 0.5 and @= 3.0 and Moo = 0.63 and a= 2 


respectively. Both the cases use 32x5 panels and 25x5x8 

field elements. These results are compared with those of 

sevi over airfoil section. All the three rectangular 
wing cases mentioned neue are in good agreement with the 


datum solutions. 
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APPENDIX A 


A.1l Evaluation of Source Influence Coefficient Bs 
Potential induced due to a unit strength of plane 


quadrilateral source element z at 7) control point 


(x,y,z) is given by 


i 


B 


d 
ij a. oe. (Aa) 


if 
j “aj 


we evaluate this coefficient in a coordinate system (§, 1,2) 


attached with the panel. In these coordinates forZ% = 0 


ff ae (A. 2) 
B,., = —_._ $$$ A. 
*J Ce ay Cae) 2 & (2) 27? 


(x,y,z) is the coordinate of control point P in the panel 


coordinate system. Change (A.2) to polar coordinates 








usimrg 
xr = ¥ R? pee and R = (62 3)° & (H= 9) 
A.3 
4 k+1 R — 
B.. = fF f [ f aS... |) ag (A.4) 
k =1 2D es (R2 422)1/2 
Performing R integration 
. 2 Oe Oi? | 
3 = op otk cp? 4 27)?72 2 a Jag (A.5) 
43 kal g 
. k 
along any side of the boundary 
R? = ae + 1? and tan J = i 
a (A. 6) 
a dl 
wo aT 
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wostitute (A.6) in (A.5) to get 
4 l 
k+1 a 
, a dl 
ae = 2 f ("+ g? 1/2 -|2 }] 
kel 1, ) + B- 
(A.7) 
ere 
2 . 2 2 
g“ = aX +2 (A.8) 
reak (A.7) in partial fractions 
| 4 k+1 2 -2 
Bay = 2 f eoeninneiane + —<= 5 +121 —ja1 
k=1 ly ce +g )1/2 (1? +a 2) (2 +g ye ae +a’) 
(A.9) 
4 Ss > : 
= E[8 t,+42°81,-8 iztz,] (A.10) 
k=1 
4 | ~ 2 
= > aI, + 2° H (A,11 
zt t J 
here 
as a 
H = a 1. - lz I. : (AT 2) 
nd 
k+1 
dl , 
ee te ———_——_. (A.13a) 
a 1, Ci + gJi/2 
Mead dl 
I= wre nen (A.13b) 
f 9. 
ne 1 (17437) (1749%)1/2 
"k 
a 
ee: oe on (A.13¢) 
3 2 -2 
l (1". + a”) 
“k 


can be easily integrated to obtain 





for k = 1, and k+1 = 2 





2 2 
Vizag + 1, 


2 
Qj +g 
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(A.14) 


I, = in ? 7 5 ; oo: 1,» 1, > Oo (A.15a) 








“4 
2 2x17 2 
= In : a Spee is dey 1G (A.15b) 
2 22372 a Gey \the - 
(1s + g*) ae 
1/2 
2 z 
[ (lf+g*) “1 Cas + g7)172 4 1] 
= in i - il > 
g 
(A.15c) 
‘. "% 1 fea? 1 IZ yf k+1 
IZl a (1 24g2) 172 lal |, 
k 
(A.16) 
1 ~1 1 ted 
3 lal laloly 
| k 
Use Eqs. (A.12), (A.16) and (A.17) to get 
= = n 5 io = 
IZllalt | (1° + g¢) a ens 
| k 
- 2 lan a 
2 Se Se 
IZtta | g* +/Z! (149°) / 1 
(A.18) 
Use the following relation (A.19) in (A.18) 
m -1 = tg * (A.19 
—— tan |m]j| = tan’ m (A.19) 


[m 
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= 1 














_ a l k4+1 
H = ae tan7? ( ) 
IZ | 2 2 2.1/2 
Q°+1ZH(1" 4# g%) k 
for KK = 1 
oe atl, apa ee ES ~1 24, 
. [Zz | [tan a) ~ tan 1 ] 
2 Cc 
1 
= t 
: a. Ula Ce @ 1. .e) 
H = oe an tan d [ Ls t 7 Z | 
' t Be 
Cy C5 + a 14 il, 
(A.20) 
where | ) | 
C5 = g? + {ZI as + g2) 172 
Ci = g? +1 ZI eS + g2) 1/2 
substitute (A.20) in (A.11) to get Bs : 
4 ! 
a (l, ec! - 1, ec!) i 
Bis = EF { a r- IZ I tan”* [ ———4+., +4 _3} : 
— t t ‘ i 
k=1_ ec’, + a 1,1, 
(A.21) 


A.2,. Evaluation of Doublet Influence Coefficient Cy. 
i sa ea cee ee a a 


Potential induced due to a unit strength of plane 
quadrilateral doublet element (25) at i a control point 


P(x, y,z) is given by © 
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das (A. 22) 


This coefficient also is evaluated in a panel coordinate 


system in which normal always points towards the +ve Z 


direction 
(CE = x)" + (tpy)~ + 2%) 
a& dan 
= 4 ie ee a (A. 23) 
| {( fx)" + (ty) + 2°} 


where (x,y,z) is the coordinate of control point FP in 


local (Panel) coordinate system and (&,71, 0) is a point 





on the panel itself. 


Change (A.23) to polar coordinates using (A.3) 





| | g R 
. | 4 k+1 
k=l BY O (R2 4 2237/2 
(A. 24) 
Performing R integration 
| 4 okt pa | as 
C., = 2. jz) 7 7. /23 % 
| (A.25) 
Using (A.6) in (A.25) 
: | a al 
ask. [ reed ich 1 , foe 
i a oe. [Z| a? 4 g*)i/2 124 32 





: a=, 4 _ = k+1 4 
1 ee Mee Gia a = 2 Ef (+H) 
| 2 j k=l 
k 
(A. 26) 
substitute (A.20) in (A.26) to get 
4 a (i | 7 
Z 7 a Cc, —- l.,¢ 
Cl..oUs sas 5 tan 1 psi i 2, (A.27) 
*J } et Kel ecole e- ty, 4 
fp “Se 


While evaluating (A.21) and (A.27) proper sign of 4, 1, 


and 1. has to be taken as shown in Pig. (A.1). From the 


Fig. (A.1) sine: and cosine of the side 1,2 w.r.t. X axis is 





Cc, = , : - = (nN, -1,)/4, (A. 28) 


og pe 
where d, =[(5, - §,)° + M, - )) ] is the length of 


the side 1,are lengths 1, and 1, are given by 


1 (§, ~ x) C, + (NW, - vy) 5S, 


Hl 


(A.29) 


+ (1, - y)S, 


Ls. (by =X) Cc, 


signed perpendicular distance a to side 1,2 is given by 


a = (§& —-x) Ss = M, - y) Cy 
| : : (A.30) 
or 


a = (5. - x) S; - (nN, - y¥) Cy 





43 


Far Field Approximation 
If the control point P(x, y,z) is too far from the 


source and doublet panel ae whose influence is sought at 


this control point P, then we can approximate the distri- 
buted singularities on panel b, by point singularity 


at its average point. 


Let t be the maximum dimension of the panel ®, 
(i.e. length of the largest diagonal) and |Rj| be the 


distance between the control point P and the average point 


of the panel vy 


R, = (x-6)F + y-q) 9 + 2 -)% 
(A631) 
where (x,y,z) is the control point P and (5 i is 
the average point of the E panel in global coordinate 
system. 
Hess an smith* have shown that if IRI /t > 4.0 we 


can use the far field approximation. In this case 


dA A | 
B = {f—- = ——— (A. 32) 
ij ee Taj ERI | 
and | 
C gees ye ek we eee ae 
ij ‘“* @n Ry * 
C. 2 = woe a 33 
3 
4 IRoI 


where A is the area of. Ds panel 


and M is the unit normal at the average point of the ay panel. 


e 
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APPENDIX B 


Evaluation of Wake Influence Coefficient Cc 
Ww 


ees 


Potential induced due to a unit strength of wake 


= th | 
element 2) at i” control point (x,y,z) is given by 


re . 
W.. 2 = Te \—=——)_ d& 
ij ws a ij = 
Y, = @ . 

= f f n ° “(— ) ag dj (B.1) 

a 4 ij 

Yo TE 

i. ee (B.2) 


~ (x— £) 2 + (y-n)? + (Zz -2) 


where (&7, 2% is a point on the wake element yin 
| j 
global coordinate system. Assuming the wake to be in Z= 0 


plane 


+ (y= 


wn | 2 - d | 1 
nN .e ee = oo ee 0 
uF Ty 8 CG? i? =e @at)4) 


eae esc l (B.3) 


f(x)? fe (y-n)? + 27177 
use (B.3) in (B.1) to get 
adgdan 


ij XY, Xpg £Geb)? + (yen)? 4 27937? 


(B.4)- 
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Performing the & integration with substitution x-§ = A 
and (y ee a ae eG 
"p Z (fx) - 


Cc = f 
Wig Ya ‘Son y) og oli GE) + (t- y)* poe ras) * 

















*DE 
-% 2 . 
rn 
y [ (-y) fay 
a 
(B.S) 
from figure Xo = Cy. + 71 tany (B.6) 
Use (B.6) in (B.5} 
Y 
b x= (C. +rn tany) 
ee | cee 6 4 ees . 
"43 Vx [(n-y)* 4 2° ] ¥ (x-C = Ntanv) 2+ (n-y) 2472 
(B.7) 
For a rectangular wing y = 0 => tanY¥Y = 0 
oy (X = C_) 

C = i é {1 + = : 3 or, 
a [(N- y)2 + 2°] ((x-C_) “+ (n- yee) 
Wij 1 2 | 

y . 2 ) 
where I, = / 2 an 
| n= + 
| | ¥.: (7 ) Zz 
1 Ney yb 
Pinot, ean ser ds (B.9) 


and 
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*b Z( x = Cc.) 
vy [ome 0? se ]ptr yoey eee 
(B.10) 


make the following substitutions in succession to integrate 








(B.10) 
H - Y =A _ and (x -C_) + 2° = 37 
ROT CA + B*) = +? 
and 
Z 
t = a ie tane 
(x - CL) (nm = Y) 
we get I = Se ea 
Z “3 Z ((n-y)? + (x-6.5 42°)" Z 
Ya 
(x - C,) (n- Y) Yp 


Ya 


of n-y 
Cc = ¢t + t 
Way | 0 z an a > f (n-¥) Die. a Eye +222 ) 


(B31) 


For a swept wing tan Y #0 and Henee ene closed form solution 
of integral (B.7) is quite complicated and consumes too much 
computer time. So we divide each ay into two parts. One 

is a semi-infinite sti emanating from the modified trailing 


edge at x = X py Contribution of this part is evaluated 


using integral (B.11). The other part consists of a flat 
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quadrilateral panel in the z = 9 Plane as shown in the 


Figure B.e2 The contribution of this part is evaluated 


using (A.27). Finally both the contributions are 


* 


summed to get the coefficient C. 
| ij 











48 


APPENDIX C 
EVALUATION OF VOLUME ‘SOURCE INFLUENCE COEFFICIENTS 
_—— A Reena than m—e nnn nesteartnneonnenteessenn 
Potential induced due to ae volume element Vs CS. 42) 


of unit source strength at its control point P(x,y,z) in 


the flow field is given by 


av a€&andad 
Pe = J = = SSIs = = Is ineiintiesinccetenainciitotea o 
Vi. ys avs ((x- §)° +(y- 1) “+ (2-8) ~) ( 
(Cel) 
Substitute x -E= A > d& = -da 
y-~-n= B => dy * -dB (Cs?) 


we will use 


fe me CP ae 4 2)1/2 
— ae « p2y2/2 
{(C.3) 
Bo o= (a2 4 c2)172 
F x (Be 4 o2yi72 
Use (C.2) in (C.1i) to get 
I= ff{ds sy } aB aA (C.4) 
(a +B 240%) tf2 
Perform the integration w.r.t. C using the substitution 


C Ste (a24p2) 1/2 tan@ to get 
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Co 











‘| Ao $B? 4 eo 4. e | 
I= Jf er eee dB] dA =/I, dA 
ae + Be 
| (c.5) 
where 
(A cle 
I, =/1n 5 QB = fln HS) dB 
| (AM +4 Be yi 
(C.6) 
integrate by parts and use (C.3) 
2 
B“C 
I, = Bln (FES +f 5 4B 
H D 
HA Be + a2 - a2 
I, = B In aim + Cf 5 aB 
HD 
Ht  @B 2 aB 
I = (SS) + Cl eee A 
1 | (a24+B74+02 1/2 (a°4+B7407) (A°+B7) 
(c.7) 
I, = Bin ae) + C 1n ( ate) + A tan7* = (€.8) 


last term in Eq. (C.7) is integrated using the following 





substitutions in succession B = 1/t, 
2 2 
1 + (a? 4.20") =.u, 


—u = C/A tang 


wa 


2 t 2 72° NS + a’p* wares, BC" 
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Using (C.8) and (c.5) 


H+C 
lT=Bysln Ban + cy in EB an +fA tan7* AH 


Bo aA (Cc. 9) 


First and second terms in (C.9) are same as I,. These are 


integrated in the same manner as 1, to get 


I = ABln (EES) + Bc 1 ) + BP tan ~1 BH 


/H+B 
= + AC ln (==) 


tA 2 -1 HC ~]1 AH 
+ BC ln (==) + c% ¢,; ee — 
———) an aa +A tan Sa «CGA 
| | (c.10) 
Let last term in (C.10) is t. 


_ -1 ~AH | | | 
ier ae tan  «~Be dA (C.11) 


integrate I, by parts 


| 4 ase 
-1 aH RC (2 A* + A*BY + A%C*) 


A“ 43° + ok (B ct tate aR tea ome 


a ee 
(ae as Be? 6251/2 
a’ (BY + Gt). 2D Be ct 


arg a's tac + B*c* 





| dA 





1 
2 4 Be + C%) 


Ren et | | 
am 5 ‘8m ‘Be - BC 


GA + 
1/2 | . 
(A . 


i | 2.2 
1 Ae (Bac) > Gee 


(at +B $e) 72 


express the third term in partial fractions 




















CENTRAL LISRaRy 
went ee KANPUR ae 
, Al cc, No. m~* AO oy 62 ao i 
2 
A -1 AH . 2 
i. = ray tan Ba 7 BC in (EEA) + = Y i rrersscasec = cA 
| a ate B +07 (A%4e7) 
Ee 
BY @A 
+ BE, 


(Ae + B? 4 0717? (22482) 


the last two terms are same as the last term of (C.7) so 


integrating in the similar manner. 


2 
AL oe AH H+A 20 as 
oO 3 «tan Ba ~ BC In (BES) ~ 07/2 tan”? HE 
- B tan? 23H 
2 AC 
(C.12) 
using Eqs. (C.12), (C.11) and (C.10) 
dit HEA Hee Be -1 BH 
I= AB ln (=~) + BC ln(=—=) + ac ln(==>) + 5- tan ae 
eat Ms a at dea 
ay 2 BC 2 AB ae 


I has to be evaluated at all the eight corners of the cuboid 


to get the value of the influence coefficient Aaj 


Ayj = 1(2,2,2) + I(t,1,2) + Tle2e1) + 1(2,1,1) = TM, 2, 2) 


2,073) = 124) = te | (C.14) 
where 2 corresponds to &, No$o 


and 1 corresponds to § 1, C4 














Far Field Approximation 


If the control point P(x, ye2) is far away from the 


volume element NG whose influence is sought at this control 


point P, then we can approximate the distributed source 


singularity in the volume ua by a point source at its 





centroid. 


Let t be the maximum dimension of element vy (1.6, 
the length of its diagonal) and Ro be the distance between 
the control point P and centroid C$. Toe gy) of element V., 


than 





t = (CE, 859% + (ny =n)? 40 =2,92972 (cr 15a) 


Ro = (Ge = 85)? = (yan)? + (2 =4,)2)°7? (crasb) 
and volume of Vj 
v= (&, -§,) - -n,) * (83-2) (C.15c) 


It has been observed by numerical experiments that for R/t > 2.5 


for field approximation holds very good. For such values of 


R/t 





(b, 26s Ss 
A 72 


ij 3 ((x-§)? + (y ~7,)° + (2-2 )°) 
(C.16) 
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APPENDIX D 


EVALUATION OF THE SURFACE INTEGRAL OVER THE SINGULAR SURFACE o 


If the point P lies on the body surface then $! = = 
ig singular there. So we exclude this point from the domain 


Dby a small hemisphere o (See Fig. D.1). We have 


Q 1 : | | 

I, = ? Sn e) uid. (D.1) 
0 1 | 

I, = He — (=) do | _ ADe2) 


If we let go become very small then 


dg = rc ag (D.3) 


where for a hemisphere 


ff ag = 2% | (D.4) 
oO. | 
over the surface gd 


2 
n 


( 


Kir 


3 fl = | 
) = = &) = _s (D.5) 


© 


using mean value theorem we evaluate the integral | 


I, = lim Sf 9 = (=) ag 
oO — 
| . 2. 
lim JSf—5 40 
i CP) mean 70 ee 


r (D. 6a) 
ie 27 Pp 
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op 4 
o 70 
* ag 
_ ea limit, 7 
~ ‘On’ mean ro a 
is (22) tim pp pag 
on’ mean Cs 
I. = O (D.6b) 


Note: (1) If the point P lies outside the surface S in 
region D then g will be a sphere and /ffdQ = 4nx. 
oC 


Proceeding in the similar fashion 
I, = -4 TP, (D. 7a) 
I, = 0 (D.7b) 


(2) If the point P lies inside the surface S in 
region D; then o = = will be continuous 
everywhere in the region D so there is no 
need to take out og. Consequently I, and 


I, both will vanish. 
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APPENDIX E 


INTRODUCTION TO COMPUTER PROGRAM 


This computer program. solves the full potential 
equation for transonic subcritical flows. The program 
consists of mainly three parts. In the first part it 
generates the surface panels and field elements, while 
in the second part it evaluates various influence 
coefficient matrices and in the final part it solves the 
full potential equation iteratively for perturbation 
velacity potential and evaluates the velocity components 


and pressure coefficient at the surface control points. 


Input Parameters 
Input of the program consists of the following 


parameters: 


ae ASPECT - RATIO 

Le TAPER — RATIO 

3. THICKNESS ~— RATIO 
4. ROOT — CHORD 

56 _ SWEEP - LE : Leading edge sweep of the wing 
6. MACH = NO : free stream Mach No. 
16 ALPHA : Angle of attack (q) | 
8. M-LINES: Number of lines of constant X/C on the surface 
9. N-LINES: Number of lines of sonstant 2Y/b on ed surface 
10. AIRFOIL-TYPE: 1 if airfoil section belongs to | 


NACA OO Series . 


s 3 for any other airfoil shape. 
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11. SPACING-CODE-X s 1 for cosine spacing along X direction 
| 2 for equal spacing along X direction 
12. Spacing-Code-Y : 1 for half cosine spacing along 
spanwise direction 
; 2 for equal spacing along spanwise 
direction 
13. Wing-type : 1 for rectangular wing 
>: 2 for swept or tapered wing 
14. Case | : 1 for nonlifting flows 


: 2 for lifting flows. 


15. Nxg : Number of field elements along X axis on the > 


root chord 


16. Nxgu : Number of field elements along X axis ahead 
of leading edge 

17. Nxgd | : Number of field elements along X axis behind 
the trailing edge 

18. Nyg s Number of field elements along y axis 

19. Nzg : Number of field elements along z amis apove 


xy plane. 


All these inputs have to be specified in one line in the input 
file. If airfoil-type is 3 then M-lines values of X and Z 


coordinates are to be specified through the input file. 


Output of the Program 


Output of the program is the coefficient of pressure 


Pp 


Cat all the surface control points and is stored in an array CP, 
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Using the input values, the program first generates the 
surface grid, surface panels, field grid and field elements. 
Coordinates of the nodes of the surface and field grid are 
stored in arrays COORD and PG ~ COORD respectively. Coordinates 
of the corner points of surface panels and field elements are 
stored in arrays PANEL and GRID. Now, the Prandtl Glauret 
transformation is applied to all the geometrical parameters 
of the surface and field grid. This is done by dividing 
the xX eseeanetes Of the arrays PANEL and GRID by a factor 
B (p= (1 SMe Ce while the y and z coordinates remain un- 
changed. By applying the Prandtl Glauret transformation to the 
poisson equation at low. subsonic Mach numbers, we geta lineal ; 
Laplace equation whose solution serves as the starting point 
for the ieeeaeion procedure. After this, the average point of 
surface panels and centroid of field element are evaluated and 
stored in arrays POINT-AV and PG respectively. 

After defining the geometry, various influence coeffici- 
ents are evaluated and stored in matrices. These coefficient 
matrices are: 

1. Vol-Surface: influence of field sources on the surface panel: 

2. Doublet-coefficient: influence of surface doublets on the 
surface panels 

3. Source-coefficient: influence of surface sources on the 
surface panels. 

4. Double-Vol : influence of surface doublets on the 
field elements 

5. Source-Vol =: infltence of surface sources on the 


field elements. 
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6. Wake-coefficient: influence of wake on the surface panels. 


7. Wake-grid * 2 influence Of wake on the field elements. 


While evaluating the source influence coefficient matrices 
we evaluate the source strength of all the panels using external 
Neumann boundary condition and store them in the array SOURCE- 
STRENGTH. After this the above mentioned matrices are assembled 
adscva tag: XS the scheme given in Section (3.3) and stored in 
Array DOUBLET-COEFFICIENT. The matrices source coefficient 
and source-strength are multiplied to get a column vector RHi1. 
The column vector RH1 is augmented to the array DUBLET- 
COEFFICIENT to get the (n x n +1) matrix which is solved for 
perturbation potential g using subroutine aNUssiohesé-eciucica 
vector is stored in the array PHI. This perturbation potential 
g is differentiated in the chordwise and spanwise directions 
using finite difference scheme to find the velocity components 
and coefficient of pressure C,, at all the surface control points 


Cp values are stored in the array CP. 


Subroutines 


The following subroutines are used by the main program: 





dis Subroutine NACA-OO-SERIES (X, Thick, 2 
For any X value and thickness ratio t/e this subroutine 


finds the Z coordinate of the NACA-OO0O-series airfoil section. 


a3 


Input Parameter: (i) X = The X coordinate of the point on 


the surface whose Z coordinate is to be evaluated. 
(ii) Thick = Thickness ratio of the airfoil section. 


Output Parameter: Z = The evaluated Z coordinate. 


2. Subroutine Unit-Vector (VECT1, VECT2, XMOD) 
It finds the unit vector VECT2 of a given vector VECT1. 


XMOD is the modulus of the vector VECT1. 


3. SUBROUTINE DOT-PRODUCT (VECT1, VECT2, DOT) 
It finds the dot produet DOT of two given veetors VECTI 


and VECT2. 


4. Subroutine Cross-Product (VECT1, VECT2, VECT3) 
It finds the cross product VECT3 of two given vectors 


VECT1 and VECT2. 


5. Subroutine Local-Transformation (ZAI, ZEETA, ETA, PO, PB, PLOCAL) 
This subroutine transforms a vector in global coordinates 


(X,Y¥,Z) to local coordinates (&.n,%) attached with the panel. 


Input Parameters: 
(i) ZAI, ZEETA, ETA = Unit vectors that constitute the right 


handed orthogonal system attached with the panel. 
(ii) PO = origin of the panel coordinate system (i.e. average 
point of the panel). | 
(iii) Pb = Vector to be transformed in panel coordinate system. 
Output Parameter: 
PLOCAL = The vector transformed into the local 


coordinate system. 
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6. Subroutine Vol (PC, ZAII, ETAA, ZEETAA, COEFF) 
This subroutine evaluates the influence coefficient of 


a vOlume source at any point either in the surface or in 


the field. 


Input Parameters: 

(i) PC : The coordinates of the control point where the 
influence is sought. 

Cae ZAIL, ETTA, ZEETAA: Lower and upper bournds of X,Y,2Z 


coordinates respectively in the cuboid. 


Output Parameters: COEFF: The evaluated influence coefficient 


Ay: 


7. Subroutine Surface Coeff (LOCAL, CONTROL-POINT, C,b) 


This subroutine evaluates influence coefficients at 

a control point either on the surface or in the field due to 

‘gource and doublet distribution on a surface panel. 

Input Parameters; 

(i) Local: Vector containing coordinates of the corner 
points of the panel in panel coordinate system. 

(ii) eeeeen estes Coordinate of the control point in the 

| panel coordinate system where influence of 


surface singularities is sought. 


Output Parameters: 





(i) C m The doublet influence coefficient Cij 


(fi) b = The source influence coefficient Baye 
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8. Function Wake (X,Y,Z, Xte, ATA) 
This function evaluates the wake influence coefficient. 


for a wing whose trailing edge sweep is zero. 


Input Parameters: 


(i) X,Y,Z2 = X,Y,Z coordinates of the control point where 


influence coefficient is to be evaluated. 


(ii) xXTe X coordinates of the trailing edge of the wing. 
(iii) ATA = Y coordinate of end points of the wake 


strip. 


9. Subroutine Gauss (C, n, n error, PHI) 


af 
The subroutine solves a system of linear equations 


using gauss elimination technique. 


Input Parameters: 
(14) Cc = The augmented matrix which is to be solved. 


(ii) mm = number of rows of the matrix C 
2 


(iv) Error = some -~ error limit to be specified for the 
| =1 
5 


(iii) n, = number of columns of the matrix C 


purpose.of pivoting (i.e. 10 


Output Parameter : 
PHI 3 The column matrix which contains the solution vectors. 
10. Subroutine FDMX (PG, VAL,DER,Ngz, Ngx, Ngy, Ng) 
This subroutine finds the first derivative along X axis 
of parameter VAL in the field using ane order accurate finite 


aifference scheme. 
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Input Parameters: 

PGi i= Vector containing coordinates of the centroids 
of field elements. 

VAL = Parameter which is to be differentiated (i.e. 
perturbation velocity potential) 

Ngz = Number of field elements along Z axis above the 


xy plane. 


Ngx = Number of field elements. along xX axes 
Ngy = Number of field elements along y axes 
Ng = Total number cf volume elements to the right of 


the XZ plane . 


Output Parameter DER : The vector containing first 


derivative of VAL along x axis in the field. 


11. Subroutine FDXX (PG, VAL, Drexz, NGz, NgX, Ngy, Ny) 

This subroutine finas the i derivative aiong X 
axis ot parameter VAL in the tield using a9 Order accurate 
finite difference scheme. input parameters ot this subroutine 
are same as those ot subroutine FDMX. The output parameter 


is DER2, the vector containing 2nd aerivative ot VAL. 


12. Subroutine FDMY (PG, VAL, DER, Ngz, Ngx, Ngy, Ng) 
This subroutine finds the 15% derivative of VAL in 
the Y direction which is stored in DER. Input and output 


parameters are same as those of subroutine FDMX. 
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13. Subroutine FDMZ (PG, VAL, DER, Ngz, Ngx, Noy, Ng) 
It finds the ie derivative of VAL in the z direction. 


Input and output parameters are again same as those of FDMX. 


NACA~OQOQO-SERIES AIRFOIL SECTION: 


Z ordinate for any X value for such airfoil section 


is given by 


Z = +5 t/e [0.2969 Vx — 0.126K = 0.3516 x? 


+ 0.2843 x3 = 0.1015 x* ] (E.1) 


Cosine Spacing: 


Here we divide a semi-circle in angular intervals of 
equal sizes and the projection of the foot of the line 
originating from the center O at angle 8 from the X axis is 


the X coordinate. (See Fig. E.1). 


N = total number of nodes, 
; Tt 
xX, = 0.5 (1 = Cos 9,); i=l, N 


Half Cosine Spacing: 
Here only one quadrant of a circle is divided in equal 
angular intervals (See Fig. E.2). 


N = total number of nodes, 
r 


2(N - 1) 


K 
i 


Cos (@,); i= i, N ; 


1. 


6 


T° 


Bo 
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Ga 


li. 


64 


REFERENCES 


Hess, J.Le. and Smith, A.M.0O., "Calculation of Nonlifting 


Potential Flow about Arbitrary Three Dimension&l Bodies", 


ES 40622, March 1962, Douglas Aircraft Co., Long Beaeh, 
California. 


Rubbert, P.E. and Saaris, G.R., "A General Three Dimen- 
sional Potential Flow Method Applied to V/STOL Aerodynamics", 
Society of Automative Engineers Journal, Vol. 77, Sept. 1969, 
Ppe 44—~51. 


Morino, L., Chen, L.T. and Suciu, E.0., "Steady and 
Oscillatory Subsonic and Supersonic Aerodynamics around 
Complex Configurations", AIAA Journal Vol. 13, No.3, 
1974. 


Tseng K. and Morino, L., “Nonlinear Green's Function Method 
for Unsteady Transonic Flows," ATIAA Progress in Astronautics 
and Aeronautics, Transonic Aerodynamics, Vol. 81, edited by 
D. Nixon, AIAA, New York, 1982, pp. 565-603. 


Johnson F.T., James, R.M., Bussoletti, J.E., Woo, A.D., 
and Young D.P., "A Transonic Rectangular Grid Embedded 


Panel Method", AIAA Paper 82-0953, 1982. 


Erickson L.L. and Strande S.M., "A Theoretical Basis for 
Extending Surface-Panelling Methods to Transonic ELOw 


 AIAA Journal Vol. 23, No. 12, pp. 1860-1867. 


Sinclair, P.M., "A Three Dimensional Field-Integral 
Method for the Calculation of Transonic Flow on 

Complex Configurations — Theory and Preliminary Results", 
Aeronautical Journal June/July, 1988, Paper No. 1482. 


Kellog, O.D., “Foundations of Potential Theory", Published 
from Springer-Verlag, Berlin 1967. 


Maskew, B., "Prediction of Subsonic Aerodynamic Char- 
acteristics: A Case for Low Order Panel Methods", 
Journal of Aircraft, February 1982. 


Johnson, F.T., “A General Panel Method for the Analysis and 


‘Design of Arbitrary Configurations in Incompressible Flows", 


NASA Contractor Report 3079. 


Ashley H. and Landahl M. (1965): "Aerodynamics of Wings 
and Bodies", Addison-Wesley. 


sO 


13. 


14. 


L506 


16. 


17. 


65 


Suzuki, S. and Washizu, K., "Calculation of Wing body 
Pressures in Incompressible Flow using Green's 
Function Method, J. Aircraft, 17, 326-331 (1980). 


Sytsma, H.S.,. Hewitt, B.L. and Rubbert, P.E., "A 
Comparison of Panel Methods for Subsonic Flow Compu~ 
tation", AGARD-AG-241. 


Salmond D.J., "A Cell-vertex Multigrid Scheme 

for Solution of the Euler Equations for Transonic 
Flow Past a Wing: Lecture Notes in Physics, Vol. 264, 
pp. 549-553. 


Sells, C.C.L. — Plane Subcritical Flow Past a Lifting 
Airfoil, Proc. Roy. Soc. A —- 308, 1968. 


Maskew, B., Rao, B.M., and Dvorak, F.A., "Prediction 

of Aerodynamie Characteristics for Wings with Extensive 
Separations", Paper No. 31 in Computation of Viscous- 
inviscid Interactions, AGARD-CP-291, Feb. 1981. 


Maskew, B., "Influence of Rotor Blade Tip Shape on Tip 
Vortex Shedding - An Unsteady Inviscid Analysis", 
Proceedings of the 36th Annual Forum of the American 
Helicopter Society, Paper 80-6, May 1980. 





FIG.2-1 The definition of the three - dimensional body 
surface and the region - D 
= 
»P 
$,D 
A 
2. 
S n 
) B | 
ny Pj Oi = 
A 
A 
Nn 
SS Spt Sywr = 
FIG. 2-2 


Nomenclature for internal dirichlet formulation 


oT 





FIG.3-1 Definition of 1 & J lines in the 
surface grid 








a | \ P(let, Je ) 
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FIG.3-6 A volume element with its corner points 
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FIG. 3.7 Definition of geometric quantities for the differentiation scheme 





FIG. 3-8 Points on the surface where forward / backward 
differencing is used 
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FIG. Al Definition of geometric quantities of a 
quadrilateral element 
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FIG. Bt Geometry of a horse-shoe vortex emerging from 
the trailing edge 





FIG. B 2. Geometry of the wake panel for the swept wing 





FIG.D1 Exclusion of the singular point P 
from the body surface 
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